scplotter to work with 10x Visium data prepared by Seurat¶

Go back to scplotter documentation: https://pwwang.github.io/scplotter/

Mouse Brain: 10x Genomics Xenium In Situ¶

See: https://satijalab.org/seurat/articles/seurat5_spatial_vignette_2#mouse-brain-10x-genomics-xenium-in-situ

In [1]:
suppressPackageStartupMessages({
    library(Seurat)
})
# Load the scplotter package
# library(scplotter)
devtools::load_all()
# devtools::load_all("../../../plotthis")

path <- "data/xenium_tiny_subset/outs"
# path <- "data/Xenium_Giotto_workshop"
# Load the Xenium data
xenium.obj <- LoadXenium(path, fov = "fov")
# remove cells with 0 counts
xenium.obj <- subset(xenium.obj, subset = nCount_Xenium > 0)

xenium.obj

ℹ Loading scplotter
Warning message:
“cells did not contain a segmentation_method column. Skipping...”
Genome matrix has multiple modalities, returning a list of matrices for this genome

Warning message:
“Feature names cannot have underscores ('_'), replacing with dashes ('-')”
Warning message:
“Feature names cannot have underscores ('_'), replacing with dashes ('-')”
Warning message:
“Feature names cannot have underscores ('_'), replacing with dashes ('-')”
Warning message:
“Feature names cannot have underscores ('_'), replacing with dashes ('-')”
Warning message:
“Feature names cannot have underscores ('_'), replacing with dashes ('-')”
Warning message:
“Feature names cannot have underscores ('_'), replacing with dashes ('-')”
Warning message:
“Not validating FOV objects”
Warning message:
“Not validating Centroids objects”
Warning message:
“Not validating Centroids objects”
Warning message:
“Not validating FOV objects”
Warning message:
“Not validating FOV objects”
Warning message:
“Not validating FOV objects”
Warning message:
“Not validating Seurat objects”

An object of class Seurat 
541 features across 36553 samples within 4 assays 
Active assay: Xenium (248 features, 0 variable features)
 1 layer present: counts
 3 other assays present: BlankCodeword, ControlCodeword, ControlProbe
 1 spatial field of view present: fov
In [2]:
options(repr.plot.width = 10, repr.plot.height = 5)

FeatureStatPlot(xenium.obj, features = c("nFeature_Xenium", "nCount_Xenium"),
    facet_scales = "free_y")

Warning message in GetAssayData.StdAssay(object = object[[assay]], layer = layer):
“data layer is not found and counts layer is used”
No description has been provided for this image
In [3]:
options(repr.plot.width = 6, repr.plot.height = 8)
# devtools::load_all()

SpatDimPlot(xenium.obj, image = "black", features = c("Gad1", "Sst", "Pvalb", "Gfap"),
    nmols = 20000, points_size = 0.1, points_palette = "Set1")
No description has been provided for this image
In [4]:
options(repr.plot.width = 8, repr.plot.height = 12)

SpatFeaturePlot(xenium.obj, layer = "counts", image = "black",
    features = c("Cux2", "Rorb", "Bcl11b", "Foxp2"),
    points_size = 0.2, points_color_name = "Expression")
No description has been provided for this image
In [6]:
options(repr.plot.width = 11, repr.plot.height = 5, future.globals.maxSize = 1024 ^ 3)

cropped.coords <- Crop(xenium.obj[["fov"]], x = c(1200, 2900), y = c(3750, 4550), coords = "plot")
xenium.obj[["zoom"]] <- cropped.coords
# visualize cropped area with cell segmentations & selected molecules
# The segmentation boundary was not loaded anyway...
# DefaultBoundary(xenium.obj[["zoom"]]) <- "segmentation"

SpatDimPlot(xenium.obj, fov = "zoom", image = "black", features = c("Gad1", "Sst", "Npy2r", "Pvalb", "Nrn1"),
    nmols = 10000, points_size = 0.1, points_palette = "Set1", shapes = TRUE)

Warning message:
“Key ‘Xenium_’ taken, using ‘zoom_’ instead”
Warning message in SpatPlot.Seurat.FOV(object, image = image, ...):
“[SpatPlot] 'shapes' is set to TRUE, meaning the same boundaries as points will be used. You may want to provide a different boundaries for shapes. Otherwise the shapes is plotted as points.”
No description has been provided for this image
In [7]:
xenium.obj <- SCTransform(xenium.obj, assay = "Xenium")
xenium.obj <- RunPCA(xenium.obj, npcs = 30, features = rownames(xenium.obj))
xenium.obj <- RunUMAP(xenium.obj, dims = 1:30)
xenium.obj <- FindNeighbors(xenium.obj, reduction = "pca", dims = 1:30)
xenium.obj <- FindClusters(xenium.obj, resolution = 0.3)

Running SCTransform on assay: Xenium

vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.

Calculating cell attributes from input UMI matrix: log_umi

Variance stabilizing transformation of count matrix of size 248 by 36553

Model formula is y ~ log_umi

Get Negative Binomial regression parameters per gene

Using 248 genes, 5000 cells

Second step: Get residuals using fitted parameters for 248 genes

Computing corrected count matrix for 248 genes

Calculating gene attributes

Wall clock passed: Time difference of 3.678288 secs

Determine variable features

Centering data matrix

Place corrected count matrix in counts slot

Set default assay to SCT

PC_ 1 
Positive:  Slc17a7, Nrn1, Epha4, Neurod6, Nwd2, Gad1, Cpne4, Rasgrf2, Rims3, Lamp5 
	   2010300C02Rik, Dkk3, Slc17a6, Pvalb, Garnl3, Cpne6, Fhod3, Plcxd2, Gad2, Tmem132d 
	   Kcnh5, Dner, Calb1, Bhlhe22, Bcl11b, Nell1, Bdnf, Rasl10a, Satb2, Arc 
Negative:  Igf2, Dcn, Fmod, Slc13a4, Fn1, Aldh1a2, Col1a1, Ly6a, Cldn5, Spp1 
	   Gfap, Nr2f2, Gjb2, Cyp1b1, Acta2, Pecam1, Adgrl4, Pdgfra, Acvrl1, Kdr 
	   Cd93, Ccn2, Cobll1, Fgd5, Sox17, Igfbp5, Carmn, Lyz2, Pglyrp1, Emcn 
PC_ 2 
Positive:  Gjc3, Opalin, Sox10, Gfap, Clmn, Vwc2l, Zfp536, Sema6a, Gpr17, Gng12 
	   Tmem163, Prox1, Adamtsl1, Dpy19l1, Cobll1, Cdh20, Arhgef28, Igfbp5, Chrm2, Sema3d 
	   Carmn, Aqp4, Fign, Pdgfra, Cspg4, Ntsr2, Lyz2, Siglech, Adamts2, Rmst 
Negative:  Slc17a7, Fn1, Igf2, Nrn1, Cldn5, Epha4, Neurod6, Ly6a, Dcn, Rasgrf2 
	   Lamp5, Fmod, Aldh1a2, Dkk3, 2010300C02Rik, Slc13a4, Car4, Nwd2, Pecam1, Gad1 
	   Col1a1, Igfbp4, Spp1, Cpne4, Rims3, Igfbp6, Acvrl1, Cpne6, Adgrl4, Calb1 
PC_ 3 
Positive:  Cldn5, Ly6a, Adgrl4, Fn1, Pecam1, Acvrl1, Kdr, Cd93, Pglyrp1, Sox17 
	   Emcn, Car4, Nostrin, Fgd5, Zfp366, Mecom, Slfn5, Paqr5, Arc, Cabp7 
	   Cobll1, Laptm5, Acsbg1, Gjc3, Siglech, Opalin, Kctd12, Ntsr2, Trem2, Sema6a 
Negative:  Slc13a4, Igf2, Dcn, Fmod, Aldh1a2, Nwd2, Col1a1, Vat1l, Calb2, Spp1 
	   Pdgfra, Gjb2, Necab2, Slc17a6, Cyp1b1, Syt6, Nr2f2, Nrp2, Dner, Slit2 
	   Col6a1, Cpne4, Spag16, Strip2, Sncg, Thsd7a, Ppp1r1b, Gucy1a1, Mapk4, Chat 
PC_ 4 
Positive:  Slc17a7, Dkk3, Cabp7, Neurod6, 2010300C02Rik, Arc, Epha4, Igfbp4, Bcl11b, Fmod 
	   Meis2, Laptm5, Dcn, Gad1, Bhlhe22, Aldh1a2, Cpne6, Rasl10a, Lamp5, Col1a1 
	   Col6a1, Igfbp6, Gfap, Satb2, Trem2, Cplx3, Gm2115, Gfra2, Garnl3, Siglech 
Negative:  Nwd2, Calb2, Slc17a6, Necab2, Syt6, Vat1l, Nrp2, Sncg, Cpne4, Cldn5 
	   Ly6a, Dner, Gucy1a1, Thsd7a, Chat, Kctd8, Tmem163, Tacr1, Tmem255a, Cacna2d2 
	   Adgrl4, Pecam1, Cntnap4, Rmst, Cbln1, Fn1, Kdr, Cd93, Inpp4b, Acvrl1 
PC_ 5 
Positive:  Gad1, Pvalb, Gad2, Rab3b, Opalin, Gjc3, Dpy19l1, Cdh13, Sox10, Garnl3 
	   Parm1, Rims3, Tmem132d, Lamp5, Neto2, Vip, Btbd11, Cntnap4, Plcxd2, Ccn2 
	   Cort, Penk, Fn1, Fhod3, Zfp536, Spp1, Col6a1, Nxph3, Sst, Dcn 
Negative:  Cabp7, Gfap, Aqp4, Laptm5, Ntsr2, Trem2, Siglech, Acsbg1, Cd53, Slc39a12 
	   2010300C02Rik, Kctd12, Cpne4, Cd300c2, Bhlhe22, Gm2115, Ikzf1, Rfx4, Necab2, Igfbp5 
	   Sipa1l3, Cpne6, Cd68, Clmn, Rmst, Nwd2, Nrp2, Orai2, Spi1, Prdm8 

Warning message:
“The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session”
12:04:15 UMAP embedding parameters a = 0.9922 b = 1.112

Found more than one class "dist" in cache; using the first, from namespace 'spam'

Also defined by ‘BiocGenerics’

12:04:15 Read 36553 rows and found 30 numeric columns

12:04:15 Using Annoy for neighbor search, n_neighbors = 30

Found more than one class "dist" in cache; using the first, from namespace 'spam'

Also defined by ‘BiocGenerics’

12:04:15 Building Annoy index with metric = cosine, n_trees = 50

0%   10   20   30   40   50   60   70   80   90   100%

[----|----|----|----|----|----|----|----|----|----|

*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
|

12:04:19 Writing NN index file to temp file /tmp/m161047/RtmpJk6gFt/file676414025e0c9

12:04:19 Searching Annoy index using 1 thread, search_k = 3000

12:04:30 Annoy recall = 100%

12:04:31 Commencing smooth kNN distance calibration using 1 thread
 with target n_neighbors = 30

12:04:33 Initializing from normalized Laplacian + noise (using RSpectra)

12:04:34 Commencing optimization for 200 epochs, with 1669008 positive edges

12:04:52 Optimization finished

Computing nearest neighbor graph

Computing SNN


Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 36553
Number of edges: 1340890

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9587
Number of communities: 28
Elapsed time: 4 seconds
In [8]:
options(repr.plot.width = 6, repr.plot.height = 8)

SpatDimPlot(xenium.obj, image = "black", points_size = 0.1)

Warning message:
“No FOV associated with assay 'SCT', using global default FOV”
No description has been provided for this image
In [10]:
options(repr.plot.width = 7, repr.plot.height = 8)

SpatFeaturePlot(xenium.obj, layer = "counts", image = "black",
    features = "Slc17a7", points_size = 0.2, points_color_name = "Slc17a7 Expression")

Warning message:
“No FOV associated with assay 'SCT', using global default FOV”
No description has been provided for this image
In [ ]:
options(repr.plot.width = 8, repr.plot.height = 7)

crop <- Crop(xenium.obj[["fov"]], x = c(600, 2100), y = c(900, 4700))
xenium.obj[["crop"]] <- crop
p1 <- SpatFeaturePlot(xenium.obj, fov = "crop", features = "Slc17a7",
    image = "black", points_size = 0.2)

# use ext argument to crop
p2 <- SpatFeaturePlot(xenium.obj, ext = c(600, 2100, 900, 4700), features = "Slc17a7",
    image = "black", points_size = 0.2)

p1 + p2

Warning message:
“Key ‘Xenium_’ taken, using ‘crop_’ instead”
Warning message in `[<-.data.frame`(`*tmp*`, , features, value = structure(list(:
“replacement element 1 has 36553 rows to replace 11872 rows”
Warning message:
“No FOV associated with assay 'SCT', using global default FOV”
No description has been provided for this image

Mini Xenium Dataset provided by Giotto vignette¶

See: https://drieslab.github.io/giotto_workshop_2024/xenium-1.html

In [22]:
path <- "data/Xenium_Giotto_workshop"
# Load the Xenium data
g <- LoadXenium(path, fov = "fov")
# remove cells with 0 counts
g <- subset(g, subset = nCount_Xenium > 0)

g

Warning message:
“cells did not contain a segmentation_method column. Skipping...”

Error in option$fn(file.path(data.dir, option$filename)) : File not found

10X data contains more than one type and is being returned as a list containing matrices of each type.

Warning message:
“Feature names cannot have underscores ('_'), replacing with dashes ('-')”
Warning message:
“Feature names cannot have underscores ('_'), replacing with dashes ('-')”
Warning message:
“Feature names cannot have underscores ('_'), replacing with dashes ('-')”
Warning message:
“Feature names cannot have underscores ('_'), replacing with dashes ('-')”
Warning message:
“Feature names cannot have underscores ('_'), replacing with dashes ('-')”
Warning message:
“Feature names cannot have underscores ('_'), replacing with dashes ('-')”
Warning message:
“Not validating FOV objects”
Warning message:
“Not validating Centroids objects”
Warning message:
“Not validating Centroids objects”
Warning message:
“Not validating FOV objects”
Warning message:
“Not validating FOV objects”
Warning message:
“Not validating FOV objects”
Warning message:
“Not validating Seurat objects”

An object of class Seurat 
541 features across 7654 samples within 4 assays 
Active assay: Xenium (377 features, 0 variable features)
 1 layer present: counts
 3 other assays present: BlankCodeword, ControlCodeword, ControlProbe
 1 spatial field of view present: fov
In [34]:
# Simple Visualization
options(repr.plot.width = 7, repr.plot.height = 6)

SpatDimPlot(
    g,
    image = "black",
    # put shapes at last
    layers = c("image", "points", "shapes"),
    features = c("ABCC11", "ACE2", "ACKR1", "ACTA2", "ACTG2", "ADAM28"),
    shapes_border_color = "cyan",
    shapes_border_size = 0.1,
    shapes_fill_by = "black",
    points_size = 0.1,
    nmols = 10000
)
No description has been provided for this image